package oneD.dynamics;

public class Lagrange {
	private int p_;
	private double[] Xi_;

	private double[] gp_; // Gauss point
	private double[] gw_; // Gauss weight

	private double tol = 1e-8;

	public Lagrange(int p) {
		this.p_ = p;
		this.initialize();
		gp_ = new double[] { -0.90617985, -0.53846931, 0.0, 0.53846931,
				0.90617985 };
		gw_ = new double[] { 0.23692689, 0.47862867, 0.56888889, 0.47862867,
				0.23692689 };
		// Gauss integration degree 9
	}

	private void initialize() {
		Xi_ = new double[p_ + 1];
		for (int i = 0; i < p_ + 1; i++)
			Xi_[i] = -1 + (2.0 / p_) * i;
	}

	public double getN(int idx, double xi) {
		double product = 1.0;
		for (int i = 0; i < p_ + 1; i++)
			product *= (i != idx) ? (xi - Xi_[i]) / (Xi_[idx] - Xi_[i]) : 1.0;
		return product;
	}

	public double getdNdxi(int idx, double xi) {
		double sum = 0.0, product = 1.0;
		boolean flag = false;
		for (int i = 0; i < p_ + 1; i++) {
			if (i != idx) {
				if (xi != Xi_[i] && !flag) {
					product *= (xi - Xi_[i]) / (Xi_[idx] - Xi_[i]);
					sum += 1 / (xi - Xi_[i]);
				} else {
					if (flag)
						product *= (xi - Xi_[i]) / (Xi_[idx] - Xi_[i]);
					else
						product *= 1 / (Xi_[idx] - Xi_[i]);
					flag = true;
				}
			}
		}
		if (flag)
			sum = 1;
		return product * sum;
	}

	public double getI(int i, int j) {
		return this.gaussIntegrate(1, i, j, -1, 1);
	}

	public double[] getI() {
		double[] I = new double[(p_ + 1) * (p_ + 1)];
		for (int i = 0; i < p_ + 1; i++)
			for (int j = 0; j < p_ + 1; j++)
				I[i * (p_ + 1) + j] = this.getI(i, j);
		return I;
	}
	
	public double[] getIt() {
		double[] It = new double[(p_ + 1) * (p_ + 1)];
		for (int i = 0; i < p_ + 1; i++)
			for (int j = 0; j < p_ + 1; j++)
				It[i * (p_ + 1) + j] = this.getIt(i, j);
		return It;
	}

	public double getIt(int i, int j) {
		return this.gaussIntegrate(2, i, j, -1, 1);
	}

	private double gaussIntegrate(int type, int i, int j, double a, double b) {
		// scale x = c+mt
		double c = (a + b) / 2;
		double m = (b - a) / 2;
		double sum = 0.0;
		double xi, f;
		for (int k = 0; k < gp_.length; k++) {
			xi = c + m * gp_[k];
			f = this.getN(i, xi)
					* (type == 1 ? this.getN(j, xi) : this.getdNdxi(j, xi));
			sum += gw_[k] * f;
		}
		return m * sum;
	}

	private double gaussIntegrateLoop(int type, int i, int j, double a, double b) {
		double epsilon = 1;
		int n = 2;
		double sum = 0.0;
		while (epsilon > tol) {
			double old_sum = sum;
			sum = 0.0;
			// divide the [a,b] to n section
			double h = (b - a) / n;
			// compute integration in each section using gauss legendre scheme
			for (int k = 0; k < n - 1; k++)
				sum += gaussIntegrate(type, i, j, a + k * h, a + (k + 1) * h);
			n++;
			epsilon = Math.abs((sum - old_sum));
		}
		return sum;
	}

	private double simpsonIntegrate(int i, int j) {
		/*
		 * compute integrate[Ni*Nj,{xi,-1,1}] numerically using Simpson rule
		 */
		int n = (int) Math.pow((p_ + 1), 5);
		double h = 2.0 / n;
		double xi = 0;
		double sum = this.getN(i, xi) * this.getN(j, xi);
		for (int k = 0; k < n; k += 2) {
			xi = -1 + h * k;
			sum += 4 * this.getN(i, xi) * this.getN(j, xi);
		}
		for (int k = 1; k < n; k += 2) {
			xi = -1 + h * k;
			sum += 2 * this.getN(i, xi) * this.getN(j, xi);
		}
		xi = 1;
		sum += this.getN(i, xi) * this.getN(j, xi);
		return sum * h / 3;
	}
}
